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The behavior of the momentum transferred to a trapped Bose-Einstein condensate by a two-photon Bragg 
pulse reflects the structure of the underlying Bogoliubov spectrum. In elongated condensates, axial phonons 
with different number of radial nodes give rise to a multibranch spectrum which can be resolved in Bragg 
spectroscopy, as shown by Steinhauer et al. [Phys. Rev. Lett. 90, 060404 (2003)]. Here we present a detailed 
theoretical analysis of this process. We calculate the momentum transferred by numerically solving the time 
dependent Gross-Pitaevskii equation. In the case of a cylindrical condensate, we compare the results with 
those obtained by linearizing the Gross-Pitaevskii equation and using a quasiparticle projection method. This 
analysis shows how the axial-phonon branches affect the momentum transfer, in agreement with our previous 
interpretation of the observed data. We also discuss the applicability of this type of spectroscopy to typical 
available condensates, as well as the role of nonlinear effects. 



PACS numbers: 03.75 Fi 

I. INTRODUCTION 

In a recent paper (l| we showed that interesting features 
emerge in the momentum transferred to an elongated Bose- 
Einstein condensate by a two-photon Bragg pulse when the 
duration of the pulse is long enough. On the basis of numer- 
ical simulations with the Gross-Pitaevskii (GP) equation, we 
interpreted those features as due to the structure of the under- 
lying Bogoliubov spectrum: the Bragg pulse can excite axial 
phonons with different number of radial nodes, each one hav- 
ing its dispersion law. Since the energy difference between 
these branches is of the order of the radial trapping frequency, 
they can be resolved only by Bragg pulses with duration com- 
parable with the radial trapping period, as observed in In 
this work, we present a more detailed theoretical analysis that 
supports our previous interpretation and gives further informa- 
tion about the role played by Bogoliubov excitations in Bragg 
spectroscopy. 

A useful insight into this problem is obtained by studying 
the case of an infinite condensate, unbound along z and har- 
monically trapped in the radial direction p. We show that the 
response of such a cylindrical condensate retains all the rel- 
evant properties needed to interpret the observed behavior of 
a finite elongated condensate. In addition to the numerical 
solution of the full time dependent GP equation, we explic- 
itly determine the time evolution of the order parameter in the 
linear (small amplitude) limit, using the quasiparticle projec- 
tion method of Refs. y|.l31- Th e cylindrical geometry allows 
us to simplify the calculations and, more important, to make 
the connection between the momentum transferred and the dy- 



namic structure factor more direct than for a finite condensate. 
By comparing the predictions of this approach to the results 
of the numerical integration of the time dependent GP equa- 
tion, one can also distinguish the different effects of linear and 
nonlinear dynamics. This analysis of cylindrical condensates 
is presented in the first sections of this paper, while section V 
is devoted to the behavior of finite elongated condensates, as 
the ones in the experiments of Refs. 01 BUI]- 

II. BOGOLIUBOV SPECTRUM OF A CYLINDRICAL 
CONDENSATE 

First we investigate the excitation spectrum of a cylindrical 
condensate. We use the mean-field Gross-Pitaevskii theory at 
zero temperature, assuming the condensate to be dilute and 
cold enough to neglect thermal and beyond mean-field contri- 
butions. The starting point is the time dependent GP equation 
for the order parameter ^(r, t) of a condensate of N bosonic 
atoms subject to an external potential V(r) |6]: 

ih W^= (-^ + V + g\m 2 )*, (1) 
at \ 2m J 

where g = AnTi 2 a/m, and a is the s-wave scattering length, 
that we assume to be positive. The order parameter is nor- 
malized according to J dr\^\ 2 — N. What is usually named 
Bogoliubov spectrum is the result of an expansion of \I/ in 
terms of a quasiparticles basis in the form Qiao]: 

*(r,t) = e-^/ h {Mr) + 

E 3 [c i u i (r) e -^* + c* U *(r)e^ i ]} ! (2) 
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where the Cj's are constants and ji is the chemical potential. 
The quasiparticle amplitudes u and v obey the following or- 
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thogonality and symmetry relations: 



dr{ 



UiU* - ViV*} = 5 



dx{uiVj — ViUj} — . (3) 



This expansion can be inserted into Eq. Q. At zero order in 
u and v, one gets the stationary GP equation: 



rv 

2m 



+ V(r) + $|Vto(r)| J ) = l*M*) 



(4) 



Among the solutions of this equation one finds the ground 
state order parameter, which can be chosen real and equal to 
the square root of the particle density 1 10]. The next order 
gives the coupled equations 

hujiUi = [H - fi + 2g^l]ui + g^lvi (5) 
-fUJiVi = [H - fi + 2gtp 2 ]vi + giplui , (6) 

where Hq = — (fi, 2 /2m)V 2 + V(r). The solutions of these 
equations provide the spectrum of the excited states in the lin- 
ear (small amplitude) regime. 

Now we treat the specific case of a cylindrical condensate, 
which is unbound along z and harmonically trapped along the 
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radial direction, p = [x 2 + y 2 ] 1 ' 2 . Let us write the trapping 
potential as V(r) = Vr(p) = (l/2)mu} 2 p 2 . The ground state 
order parameter is uniform along z, while the stationary GP 
equation for the radial part becomes 



h'v 2 



4>o(p) = p.(/)o{p) (7) 



where N/ L is the number of bosons per unit length and the 
function (f>o is chosen to be real and subject to the normaliza- 
tion condition 2ir J^°dp pcf>Q = 1. In this geometry, the quasi- 
particle amplitudes u and v can also be factorized. We are 
actually interested in axially symmetric states (azimuthal an- 
gular momentum equal to zero), since they are the only ones 
excited by a Bragg process in which the momentum is im- 
parted along the z-axis. These states, characterized by the 
axial wavevector, k, and the number of radial nodes, n, can be 
written in terms of new functions u n: k{p) and u n> fc(p) in the 
form 

(u,w)„ )fc (p,«) - L-^ 2 e lkz (u,v) ntk (p) . (8) 
The Bogoliubov equations 0-© then become 



huj n , k u n , k (p) 



2m 



2 m 



h 2 k 2 

2m 

h 2 k 2 

2m 



V T (p) 
V T (p) 



~ ? -4>1{P) | V nt k(p) + ^-^o(p) u ",fc(p) 



L 



L 



(9) 
(10) 



Numerical solutions of the Bogoliubov equations for cylin- 
drical condensates have already been obtained in the context 
of theoretical studies of the Landau critical velocity 11111, the 
Landau damping of low energy collective oscillations 1131 . 
and the behavior of solitary waves 1 13]. Here we solve equa- 
tions (I9l- (I10> in order to use the functions u and v as a basis 
for a time dependent calculation, as explained in the next sec- 
tion. 

It is worth noticing that, expressing all lengths and energies 
in units of the harmonic oscillator length, a p = [^./(mWp)] 1 / 2 , 
and energy, Tilu p , respectively, the solutions of GP and Bogoli- 
ubov equations for the cylindrical condensate scale with a sin- 
gle parameter, aN/L. The same parameter can be expressed 
in terms of the chemical potential in the Thomas-Fermi limit, 
Mtf- One has in fact, rj = ptf/(^ p ) — [iaN/L] 1 ^ 2 
llU. UM . We perform calculations for different values of 77, 
which are representative of available elongated condensates. 
The lowest value that we consider is 77 = 9.4, which simu- 
lates the condensate of Davidson and co-workers Cj H|, the 
largest value is 77 = 70, for the condensate of Ref. 11311 . We 
also use an intermediate value, rj — 26.5, which corresponds 
to the condensate of the calculations in Refs. fl^.fl7ll . 



In Figs, la and lb we show the frequency w„ fc of the Bo- 
goliubov branches as a function of k, for 77 = 9.4 and 26.5, 
respectively. The lowest branch starts linearly at low k and 
then becomes quadratic for large k, as expected for the text- 
book Bogoliubov spectrum of a uniform condensate. The 
transition between the two regimes occurs at k of the order 
of = Jr\, where £ is the healing length. The slope of this 
branch has been the object of the investigation of Ref. fiUl . 
The second branch starts at lo 2 ,q = 2ui p , where it reduces to 
the transverse monopole (breathing) mode; its frequency is 
model independent as a consequence of a scaling property of 
GP equation in two dimensions 1 18J. For k <C the dis- 
persion law of the lowest branches can be obtained, with great 
accuracy, also within a hydrodynamic approach 1 1^1 . 



Our purpose is now to explore the effects of the discrete 
multibranch spectrum of Figs. 1 a and lb in the response of the 
condensate to a Bragg pulse. 
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FIG. 1: Spectrum of axially symmetric Bogoliubov excitations of 
cylindrical condensates with r\ = 9.4 (a) and r\ = 26.5 (b). The 
frequency Lu n ,k, in units of the radial trapping frequency ui p , is plot- 
ted as a function of the axial wavevector k, in units of aT , with 
a p = [h/ '(mujp)] 1 / 2 . The number of radial nodes is n = 0, 1, 2, . . . , 
starting from the lowest branch. The shaded area corresponds to the 
region that gives finite contributions to the the dynamic structure fac- 
tor S(q, lu) in local density approximation (see section HVl . 



III. TIME DEPENDENT QUASIPARTICLE AMPLITUDES 
AND MOMENTUM TRANSFERRED 



Let us suppose that a Bragg pulse is switched on at a time 
t = 0. In the experiments, this is obtained by shining the 
condensate with two laser beams with wavevectors ki and 
k2 and frequencies io\ and u>2- The action of the stimulated 
light-scattering from atoms of the condensate is equivalent to 
a moving optical potential with wavevector q = ki — k2 and 
frequency u> — lo\ — lo-i- Within the GP theory, this is ac- 
counted for by adding an extra term to the external potential, 
that becomes 1 20] 



where 9(t) is the Heaviside step function and we have chosen 
q along z. We calculate the response of the condensate in 
two ways: i) by using the quasiparticle projection method of 
Refs. |2j|3||; ii) by numerically solving the full time dependent 
GP equation Q for given values of 77, Vb, q and ui. The latter 
case will be discussed in the next section. Here we briefly 
sketch the quasiparticle projection method and its application 
to infinite cylindrical condensates. 

The starting point is a generalization of the Bogoliubov ex- 
pansion 10 that allows the population of quasiparticle states 
to vary with time 1 2]. The simplest way is to expand the order 
parameter in terms of a quasiparticle basis as in (|2j, but taking 
the coefficients Cj to be time dependent: 



*(r,t) 



-ifxt/h 



Oo(r) 



]} -(12) 



The functions tpo, uj and Vj are the solutions of Eqs.@-(|6j 
and are assumed to be static. Let us suppose that at t = the 
condensate is in the ground state and thus all Cj(0) are zero. 
For t > the Bragg potential starts populating the quasipar- 
ticle states. The dynamics is entirely fixed by the evolution 
of the coefficients Cj(t), that can be obtained by inserting the 
expansion H2\ into the GP equation Q, with the external po- 
tential given in flit . One finds H 

3 

= ffl^ e-*"*/* cos(qz - umtr, t) . (13) 
in 

In order to get an expression for the c's, one can multiply by 
u* and v* the last equation and its conjugate, respectively, 
and integrate in dr. By using the orthogonality and symmetry 
relations fJJ one gets 



6{tW_B 
in 



x cos(qz — tot) 



(14) 



Finally, assuming that the quasiparticle occupation is always 
much smaller than N, one can replace ^> with the ground state 
order parameter and integrate in time. The results is fJI 



J 



1 giujit 



dr (u* + v*)tp cos(qz - tot') 



(15) 

which is expected to be valid in the limit of small perturba- 
tions. 

In the case of a cylindrical condensate, we can use the fac- 
torization (|8} to rewrite the previous expression in the form 



dPP[u* n ,k(P) 



inL 



V(r, t) = V T (p) + d(t)V B cos(qz - tot) , 



(11) 



+<fc(p)]<Mp) / dze ikz cos(qz - LUt') . (16) 

J —00 

The last integral simply selects out the quasiparticles with k — 
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±q, so that 

Cn,k{t) 

where 



NV B W n . k 



2h 



OJn.k — W 



, (17) 



W n , k = 2vr y dp pK :k (p) + v* tk (p)]Mp) ■ (18) 

This means that, once the ground state and the Bogoliubov 
spectrum are known at t = 0, the order parameter at time 
t > can be calculated in a rather simple way. 

An important quantity that can easily be calculated within 
this scheme is the total momentum imparted to the conden- 
sate, which can be defined by integrating the current density 
associated with the order parameter: 



r d 
- / dr **(r,t) — *(r,i)+c.c. . (19) 

i J dz 

By inserting the expansion (11 21 and using the orthogonality 
and symmetry relations (0, one gets the simple expression 



P z (t) =5>fc|c„, fc (;)| 2 



(20) 



71 , k 



Now, let us take c n ^(t) from Eq. d!7i and use the fact that the 
eigenfrequencies and eigenfunctions of the Bogoliubov equa- 
tions d9l-(fl"0l are even functions of k. One gets the final result 



sin [(a;. 



7i, q 



#iy 



(w ntq - Lj)t/2 J 



ain[(uJn,q + U))t/2] 
(bJ n . q + uj)t/2 



(21) 



In the next section we will show typical results obtained with 
this expression, but some relevant features are already evi- 
dent. Notice first that, for positive u> and a given q, the mo- 
mentum transferred d2U is clearly resonant at the frequencies 
w = Wm,ij. In other words, a peak occurs in P z whenever 
a vertical line at k = q crosses the Bogoliubov branches of 
Figs, la and lb. The separation between these branches is 
roughly 1-2 ui p , while the width of each peak goes like 2n/t. 
Thus, in order to resolve the different peaks, one has to wait 
at least a time of the order of trapping period T p = 2w/u} p . 
For t in this range, the second term in the square bracket gives 
certainly a negligible contribution. Finally, the contribution 
of each axial branch to the total momentum depends on the 
quantities W n q , defined in Eq. Jl 81 . Typically, for a given q, 
the lowest branches have a greater weight in the summation 
( 12 1> . since the corresponding quasiparticle amplitudes have a 
greater overlap with the ground state. 

This behavior of P z reflects its direct connection with the 
dynamic structure factor. At zero temperature, the latter is 
defined as 



S(q,c)=^|(z|p q |0)| 2 <5(c-^) 



(22) 



where p q is the density fluctuation operator, |0) is the ground 
state, \i) are the excited states, and uji — (Ei — En)/h. In 
terms of the quasiparticle amplitudes one has I2lt l22ll 



i|/5 q |Q)= /drK(r)+<(r)]e^ r Vo(r) 



(23) 



If q is taken along z, and the system is translational invariant 
in the same direction, one can use Eqs. (|8} and dl 8t i to rewrite 

S(q,uj) as 

S(q,Lu) = Y,N\W n J 2 6{u;-Lu n , q ). (24) 

n 

Comparing this result with Eq. J21I one finds l3l fT6ll 

lim P z (t) = l!?¥t[S( qt u) - S(-q,-u)] . (25) 
t— too zn 

As discussed in Ref. |3], these two quantities are connected 
in such a simple way only for a cylindrical condensate. If 
the condensate is trapped also along z, the relation between 
S(q,u>) and P z is less direct, involving a further integration 
on the duration time of the Bragg pulse. However, if the ax- 
ial trapping frequency is significantly smaller than the radial 
one, there exists a wide range of time where the connection 
between the momentum transferred and the dynamic struc- 
ture factor is roughly the same as for an infinite cylinder. In 
that range, the Bogoliubov branches, which correspond to 8 
peaks in S(q, ui), are also observable as distinguishable peaks 
in P z (t), as we will see in section V. 



IV. BOGOLIUBOV BRANCHES AND MULTIPEAK 
SPECTRUM 

In Figs. 2a and 2b we show the momentum transferred to 
two cylindrical condensates with rj — 9.4 and 26.5 by a Bragg 
pulse with wavevector q = 2.3a" 1 (a) and q = 4 a^ 1 (b). 
Recalling that = Jfj, the values of q are chosen in order 
to have roughly the same value of q£ in the two cases. As solid 
lines we plot the quantity P z (t) / (NV^hq) as a function of u> 
and for three different durations times: t = 2, 6 and 10 w" 1 . 
According to Eq. (12 1> . this quantity is independent of Vb- 

For a pulse shorter than the trapping period, the response is 
distributed over a single broad peak. Short pulses were used in 
the first experiments of Refs. 0,|5|]. In this regime, accurate 
predictions can be found by using the local density approx- 
imation (LDA) fl6l l22ll . i.e., assuming the system as locally 
uniform. If the excitation spectrum is the one of a uniform gas, 
but with a dispersion law fixed by the local density n(r), the 
dynamic structure factor S(q, lj) turns out to be non zero only 
within the frequency range, uj r < uj < uj r [l + 2\xj (hwr)] 1 ^ 2 , 
where Tiuj r — h 2 q 2 / (2m) is the free recoil energy |22J. This 
range corresponds to the shaded area in Figs, la and lb. The 
LDA approximation for S(q, uS) can be used to estimate the 
momentum transferred P z (t) as in |16]. The corresponding 
curves are shown as dashed lines in Figs. 2a and 2b. The 
agreement between dashed and the solid lines for lowest curve 
at t = 2u)~ x is reasonable; for shorter times we reproduce the 
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FIG. 2: Momentum transferred to cylindrical condensates with 
r\ — 9.4 (a) and r\ — 26.5 (b) by a Bragg pulse of wavevector 
q — 2.3ap 1 (a) and q = 4 dp 1 (b), and frequency u. The quan- 
tity P z /(NVghq), in units of (?kj p )~ 2 , is plotted vs. u, in units 
of Lu p . Solid lines correspond to the results obtained from Eq. <2H 
at three different times, t — 2 (bottom curve), 6 (mid) and 10 (top) 
u7 . The dashed lines are the corresponding predictions of the local 
density approximation. The triangles on the horizontal axis are the 
frequencies of the lowest Bogoliubov states that contribute to P z . 



results of Ref. fTfjl . and the agreement with LDA becomes 
better and better. 

For t of the order of T p — liruip 1 new structures begin to 
appear, strongly deviating from LDA. In particular, one finds 
a multipeak spectrum that reflects the existence of Bogoliubov 
branches, i.e, the same plotted in Fig. 1 . The peak at the lowest 
frequency is the collective axial mode with no nodes in the 
radial direction (n = 0). The n = 1 and n = 2 branches are 
also visible for the condensate with 77 = 9.4 in (a), and even 
more for 7/ = 26.5 in (b). The corresponding frequencies, 
u) n>q , are indicated by triangles on the horizontal axis. 

From Eq. (12 1 1 one sees that the height of each peak, for 
t ^s> T p , becomes proportional to t 2 \W nyq \ 2 . In Figs. 3 we 
plot iVKn^l 2 as a function of q for the same condensates of 
Fig. 2. This figure tells us how many branches give signifi- 



FIG. 3: Quantities IW^.q] 2 , defined in Eq, <18> . for the two conden- 
sates in Fig. 2. In the linear response regime, these functions provide 
the height of the <5-peaks in the dynamic structure factor 1241 . as well 
as of the peaks in the momentum transferred J2H . The values of q 
used in Fig. 2 are shown as vertical dotted lines. 



cant contributions to P z and S(q, u) at a given q. One notices 
that the contribution of each branch is sizeable when the cor- 
responding u) n>q is, roughly speaking, within the shaded area 
in Fig.l. The first branch is the dominant one at small q. Its 
weight |Wo, 9 | 2 also vanishes at q — > 0, in agreement with 
the limiting behavior of long wavelength phonons, for which 
uo :q ~ ~vo,q 9 so that the integral dl 8t i vanishes. 
One also notices that, for a given q, the number of branches 
that contribute to the summation i2l\ increases with r\. This 
effect is rather dramatic for the condensates used by Ketterle 
and co-workers [4], where r\ ^> 10 (see Fig. 4). In this case, 
even a small extra broadening in the experimental detection of 
P z , might significantly reduce the visibility of the Bogoliubov 
branches. 
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FIG. 4: Same as in Fig. 2, but with 77 = 70, q = 6 a~ x and t = 
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The momentum transferred can also be calculated by direct 
integration of the time-dependent Gross-Pitaesvkii equation 
([D, using the definition (1191 for P z . In the linear response 
regime the results must concide with the ones of the quasi- 
particle projection method so far discussed. An example is 
shown in Fig. 5. The solid line corresponds to the momen- 
tum calculated with Eq. (12 It at t = 6oj~ \ as in Fig. 2a. The 
empty squares are the results of GP simulations for a weak 
Bragg pulse (Vb — 10 -3 hw p ). As expected, the two results 
are in good agreement, except for small differences that are 
compatible with the accuracy of our GP simulations. We also 
checked that, for such small values of Vb, the quantity P z (t) 
scales with Vg as expected for the linear response. 

A major advantage of using both approaches, is that one can 
now increase Vb, beyond the linear regime, and look at the 
differences between the linear response, given by Eq. J2 1 b . 
and the numerical GP calculation, which is valid also in the 
nonlinear regime l3l l23ll . A typical result is shown in Fig. 5, 
where we compare the small Vb limit (solid line and empty 
squares) with the results of GP simulations for Vb = 0.5 hu> p 
(solid circles). With such a pulse intensity, the fraction of con- 
densate atoms which are excited is roughly 25% at resonance. 
This value is slightly larger than that of typical experiments, 
where a fraction of the order of 10-20% is required for the 
optical detection of the excited atoms. The figure shows that, 
even for such highly excited condensates, the main features 
of the multipeak response, associated with the underlying Bo- 
goliubov spectrum, are still visible. The most significant ef- 
fects of nonlinearity, originating from the mean-field interac- 
tion in the GP equation, are that i) P z does not scale with Vg, 

ii) the peak frequencies are slightly shifted downwards, and 

iii) additional structures appear superimposed to the original 
shape of the peaks. This behavior is a consequence of typical 
nonlinear processes like mode-mixing and higher harmonic 
generation. Similar effects were already predicted for spheri- 
cal 01 and one-dimensional 1 23] condensates. 

V. BRAGG SCATTERING FROM A PROLATE 
ELONGATED CONDENSATE 

Here we show that the main features found for the response 
of a cylindrical condensate to a Bragg pulse remain observable 
also in a finite elongated condensate. 

Let us take the trapping potential in the form 

V t (t) = {l/2)muo%p 2 + X 2 z 2 ) (26) 

where A = uo z /uj p . If A < 1 the condensate at equilibrium is 
a prolate ellipsoid. The ground state at t = can be found as 
the stationary solution of Eq. Q. Then, the time dependent 
GP equation can be solved at t > to simulate the Bragg 
process. 

In Figs. 6a and 6b we show typical results obtained for 
P z in the linear response (small Vb) regime. We simulate 
the Bragg scattering from two different condensates, whose 
chemical potential, in units of hui p , is equal to 9.4 and 26.5, 
as for the cylindrical condensates of Fig. 1 and 2. These two 
values correspond to the condensate used in the experiments 
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FIG. 5: Momentum transferred to a cylindrical condensate with r\ — 
9.4 by a Bragg pulse with q = 2.3 a" 1 . The quantity P z /(NV£Kq), 
in units of (hio p )~ 2 , is plotted vs. u), in units of uj p . The solid 
line corresponds to the result obtained from Eq. 12 li at t = 6 LOp 1 
(same as in Fig. 2a). Empty squares and solid circles are the results 
of the numerical integration of the time dependent GP equation Q 
for a Bragg pulse of intensity Vb = 10 -3 hu> p and Vb = 0.5 hto p , 
respectively. Dashed and dot-dashed lines are a guide to the eye. 



of Ref. UJ, made of A = 10 5 atoms of 87 Rb in a trap with 
uj± = 27r(220Hz) and A = 0.114, and to the one used in the 
calculations of Ref s. |16, 17], with A = 0.125, respectively. 
With this choice, the condensates in Fig. 2a (2b) and 6a (6b) 
have also the same density profile in the radial direction for 
2 = 0. As one can see, the momentum transferred to the 
trapped condensates of Fig. 6 has the same multipeak struc- 
ture of the corresponding cylindrical condensates of Fig. 2. 

The response of an ellipsoidal (axially confined) conden- 
sate is expected to differ from the one of a cylindrical (axially 
unbound) condensate for two main reasons. First, their Bo- 
goliubov spectra are different: the excitations of a confined 
condensate are discretized also along z and the quasiparticle 
amplitudes u and v have a nontrivial dependence on p and z. 
Second, the response for long pulses is affected by the pres- 
ence of the axial confinement, which can produce a reflection 
of quasiparticles at the boundaries and a center-of-mass mo- 
tion of the condensate in the trap. 

The first effect, however, is small if g ^> L~\ where L 
is the axial length of the condensate. In this case, the Bragg 
pulse excites quasiparticles that behave like travelling waves 
along z, having wavelength much shorter than the condensate 
size. One can safely introduce a continuous wavevector k to 
identify such excitations as in the multibranch spectrum of 
the cylindrical condensate. The axial size of the two conden- 
sates in Figs. 6a and 6b is approximately L = 75 and 120a p , 
respectively, so that the Bogoliubov modes can be safely rep- 
resented by continuous branches down to about k = 0.5 a" 1 . 
Moreover, for k 3> L , the frequencies u> n ,k for the elon- 
gated ellipsoidal condensate are expected to be close to the 
ones of the corresponding cylindrical condensate, with a sim- 
ilar multibranch structure. In order to check that the peaks 



7 




2 4 6 8 10 




5 10 15 20 25 



to 



FIG. 6: Momentum transferred to two different ellipsoidal con- 
densates by a Bragg pulse. In part (a) the condensate has n — 
9Ah(j p and A = 0.114; the Bragg pulse has q = 2.'Sa p ~ 1 and 
Vb = 0.02hujp. In part (b) the condensate has fi = 26.5Tiuo p and 
A = 0.125; the Bragg pulse has q = 4ap 1 and Vb = 0.05huj p . 
The points are obtained by numerically integrating the GP equation 
Q for different lu (in units of ui p ). The duration time is t = 2 and 
6 LUp 1 for the lower and upper curves, respectively. The lines are a 
guide to the eye. The solid circles on the horizontal axis in (a) are 
the frequencies of the lowest Bogoliubov branches for that q. 



of P z are still located at the frequencies of the Bogoliubov 
branches, we make a further simulation for the condensate in 
Fig. 6a. We excite the condensate by acting with a Bragg pulse 
for a given duration t and then we let it to freely evolve in the 
trap. We perform a Fourier analysis of the induced density 
variations in the axial direction. As expected, we find that, 
for each q and ui, the density oscillates as a superposition of 
waves of different frequencies w n The latter are extracted 
from the Fourier spectrum and reported in Fig. 6a as solid cir- 
cles. As one can see, the momentum is resonantly transferred 
at those frequencies. 

The second effect is also small if the duration time of the 
Bragg pulse is significantly less than the axial trapping pe- 
riod T z . If A is of the order of 10 , or less, there exists a 
sufficiently wide interval of time, between T p and T z , where 
the multipeak structures of P z can be resolved without being 
significantly affected by the reflection of quasiparticles at the 
boundaries and by the center-of-mass motion. 

Finally, the experimental observation of the excited atoms, 
requires rather high Bragg intensities. As already shown in 
the previous section, nonlinear processes are not expected to 
hinder the visibility of the multibranch Bogoliubov spectrum. 
Conversely, the combined use of both the quasiparticle projec- 
tion method and the time dependent GP equation might allows 
one to extract from the observed P z information about nonlin- 
ear effects in the condensate dynamics. 

In Fig. 7 the results of our GP simulations are compared 
with experimental measurements of P z . Similar examples 
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FIG. 7: Momentum transferred to the condensates of Ref. Q]. 
Empty symbols are obtained by numerically integrating the GP equa- 
tion Q for different values of u> (in units of u> p ). Points with er- 
ror bars are experimental data Q|. The lines are a guide to the 
eye. The condensate is the same of Fig. 6a. The Bragg pulse has 
wavevector q = 2.3a" 1 , duration t = 6 ms= 8.29ojT and intensity 
Vb = 0.2hu p . 



were already reported in Ref. |1], where the major sources 
of noise and broadening of the observed data were also dis- 
cussed. The main message of Ref. 1 1] was that the measured 
response to long Bragg pulses exhibits clear evidences of the 
underlying Bogoliubov spectrum. In the present paper we 
have provided a deeper theoretical analysis, which supports 
this view. 
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NUMERICAL PROCEDURES 

• i) Stationary GP equation 

The ground state configuration at t = is obtained by 
solving the stationary GP equation @ with the external 
potential given in i26\ . The case A = corresponds 
to the cylindrical condensate of sections III and IV: the 
order parameter is uniform along z and the problem is 
reduced to the one-dimensional equation @. The solu- 
tion of this equation is equivalent to the minimization 
of a GP energy functional |6]. We find the minimum by 
mapping the function (f>o on a grid of points and prop- 
agating it in imaginary time with an explicit first or- 
der algorithm, starting from a trial configuration, as de- 
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scribed in Ref. |24]. The case A ^ corresponds to an 
ellipsoidal condensate, as in section V. We use the same 
minimization algorithm. Due to the axial symmetry of 
Vt, the problem is two-dimensional. The order param- 
eter is mapped on a N p x N z grid (typically, 64 x 1024 
points). 

• ii) Time dependent GP equation 

Let us define an effective Hamiltonian H by rewriting 
Eq. ([0 in the form ihdt^ = H^>. We solve this equa- 
tion by propagating the order parameter in real time, 
starting from the ground state configuration. At each 
time step, At, the evolution is determined by the im- 
plicit algorithm 



iAt \ ( iAt , 

1 + — H\ # n+1 = f 1 - —H ) *„ (27) 



where H and are calculated at the n-th iteration. 
For A ^ 0, the propagation is splitted into the ax- 
ial and radial parts. The former is obtained by us- 



ing a fast Fourier transform algorithm to treat the ki- 
netic energy term, while the latter is performed with a 
Crank-Nicholson differencing method, as in 11611 . For 
the cylindrical case (A = 0), we first write the or- 
der parameter as $ (p, z, t) = ^2 ip v (p,t)exp(ivqz), 
where q is the wavevector of the Bragg potential, and 
then we propagate in time the radial functions <f>„, for 
v = 0, ±1, ±2, .... Only a few values of v give signif- 
icant contributions to the sum, their number depending 
on the intensity of the pulse. 

• iii) Bogoliubov equations 

In order to solve the eigenvalue problem (I9l-( I10> for a 
cylindrical condensate, we project the unknown radial 
functions u(p) and v(p) on a basis set of Bessel func- 
tions. The problem is thus reduced to a matrix diagonal- 
ization. We use Bessel functions with up to 100 nodes 
in the computational box. The radius of the box is typi- 
cally two or three times the size of the condensate. 
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